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ABSTRACT 

This paper derives the three-dimensional lambda-formulation equations for 
a general orthogonal curvilinear coordinate system and provides various block- 
explicit and block-implicit methods for solving them, numerically. Three 
model problems, characterized by subsonic, supersonic and transonic flow 
conditions, are used to assess the reliability and compare the efficiency of 
the proposed methods. 
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INTRODUCTION 


Among the many theoretical models employed in the numerical simulation of 
compressible inviscid flows the so-called lambda-formulation has received 
considerable interest (see, e.g., [1-8]): the time-dependent Euler equations 

are recast into compatibility conditions of bicharacteristic variables along 
the corresponding hi characteristic lines and discretized using windward 
^■^^^ erences 9 order to account for the direction of wave propagation 

phenomena, correctly. Such an approach has many nice properties: it provides 

very accurate numerical results, even with rather coarse meshes (see, e.g., 
[2], [3], [6]); it requires only the physical boundary conditions, so that 

there is no need for any additional numerical boundary treatments, which are 
frequently the cause of numerical instability [9]; it handles in a most 
automatic and physical ly-sound way mixed supersonic-subsonic flow fields; and 
finally, it has a well documented, although controversial, capability of 
capturing shocks without any additional dissipation [2—6]. For these reasons, 
in spite of the fact that the captured shocks" are only isentropic 
approximations to correct weak solutions of the Euler equations and do not 
correctly move within the flow field - unless properly fitted [4], the lambda- 
formulation is considered to be a very useful and reliable tool for predicting 
compressible flow fields and, therefore, very worthy of further studies and 
improvements; and in fact, in the last two years, for the cases of quasi-one 
dimensional and two dimensional flows, the development of various kinds of 
implicit integration schemes [5-8] has removed the only manor limitation of 
previous lambda methods, namely, the CFL stability restriction associated with 
their explicit integration procedures. 

It now appears very timely and worthwhile to develop efficient numerical 
methods, based on the lambda-formulation, for three-dimensional flows, as done 
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in the present paper: the "most appropriate" three-dimensional lambda- 

formulation equations are first derived for the case of a general orthogonal 
curvilinear coordinate system; the governing equations are then discretized 
and linearized in time using a delta approach and various block-explicit as 
well as block-implicit numerical techniques are proposed to solve the 
resulting discrete equations approximately at every time step; all of the 
proposed methods are finally applied to solve three model problems, 
characterized by subsonic, supersonic and transonic flow conditions, 

respectively, in order to assess their reliability and efficiency. 


THREE-DIMENSIONAL LAMBDA-FORMDLATION EQUATIONS 

The nondimensional continuity and momentum (Euler) equations for the 
homentropic flow of a perfect gas are given in vector form as [3,5] 

6 (a + V • Va) + a V • V = 0 (1) 

V t + (V • V)V + 6 a Va = 0, (2) 

where a is the speed of sound, is the velocity vector, is the gradient 

operator, t is the time, subscripts indicate partial derivatives and 
6 = 2/(y-l), y being the specific heats ratio. 

In a general orthogonal curvilinear coordinate system we have: 


1 = V 1 ®1 + v 2 *2 + v 3 


(3) 


-1 8 


-2 8 


-3 8 


v= 8 8 8 _ -- - + — 

- -1 8s 1 —2 8s 2 -3 8s 3 h x 8q 1 h» 2 8q 2 h 3 ^ 


( 4 ) 
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- • ^ % (h 2 h 3 V l ) + l^ ( h l h 3 v 2) + ( b l b 2 v 3H- (51 


where = 1,2,3) are the unit vectors of the (right-handed) orthogonal 

curvilinear coordinate system, q.. and h^ are the corresponding coordinates 

and scale factors, ds^ = h^ dq^ are the elementary arc lengths along the 
coordinate lines, see [10], and v^ are the components of V. Equations (1) 
and (2) can be written in the general orthogonal curvilinear coordinate system 

by means of eqs. (3-5) and some lengthy but straightforward algebra, the only 

difficulty being the evaluation of the derivatives of the unit vectors 
£ i with respect to the coordinates q^. These expressions are therefore 
given here as: 



e, 3h. 
h i 8q i 



e. 3h. 

—j i 


h .i 3q .i 


e, 3h . 
-k i 


\ 3q k 


(i * i) 

( i *j * k). 


(6) 

(7) 


The six lambda compatibility equations can now be obtained by summing and 
subtracting from the continuity eqn. (1) each component of the momentum eqn. 
(2), to give: 


3C 

+ v 9C + 

3C 

3s^ 

2 3s 2 

V 3 3s 3 - 

3D 

+ v 3D + 

3D 

9s l 

+ V 2 3i^ + 

V 3 3s 3 


1 P 1 ' 1 . 


E t + V 1 + ( v 2 +a ) H- + v 3 If- - - c 2 + S 2 - r 2 


F t + V 1 H; + <V a) + v 3 Hr ' 5 - “2 + B 2 + Y 2 


(8a) 


(8b) 

(8c) 

(8d) 
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_ i 3G . 3G , f , v 3G 
G t + v l ^ + v 2 ^ + ( '’3 4a) ^ 


= -c - 


a 3 + B 3 Y 3 


(8e) 


„ , 3H 3H , . 3H 

H t + V 1 3^ + V 2 3i^ + (v 3’ a) 3^ 


= C - 


a 3 + B 3 + y 3 


(8f) 


where C, D, E, F, G and H are the six bicharacteristic variables given as 


C = v^ + 6a; 


P = v x - 6a; 


E = V 2 + 6a 


F = v 2 - 6a; 


G = v^ + 6a; 


H = v 3 - 6a 


A = 


h l h 2 h 3 


K k <h 2 h 3 ) + v 2 k: (h i h 3 ) + v 3 k 'Wl 


v, v 0 3b, v 0 3h, 

1(2 1 _J3 

“1 = h l S aq 2 h 3 aq 3 J 


v 2 2 3h 2 v 3 2 3h 3 

B 1 = h^“h^ 3^ + h^h^ 3^ 


3 v, 


3 v. 


1 UV 0 V \ ~ 

f 1 2 1 

Y 1 " a ^h 2 3q 2 h 3 3 q 3 


(9a,b,c) 

(9d,e,f) 

( 10 ) 

(Ha) 

(11b) 

(He) 


and a 3 , B2’***’^3 have very similar expressions, which can be obtained by 
appropriate subscripts rotation and are thus omitted for the sake of 
conciseness . 

Equations (8) are the compatibility conditions of the bicharacteristic 
variables along their bicharacteristic lines (in the four-dimensional q 3 , 
q 2> 93* t space) obtained by the intersection of the bicharacteristic conoid 
(associated with a point P) with the q^t, q 2 ~t and q 3 ~t surfaces passing 
through its vertex P (the left-hand sides of eqns. (8) clearly being total 
derivatives along such lines). Therefore, they could be integrated by means 
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of any numerical method using windward differences according to the direction 
of wave propagation along the bicharacteristic lines, thus providing a three- 
dimensional lambda scheme. However, like in the two-dimensional case [6], 
there are two major difficulties associated with solving eqns. (8) 
numerically. First, the six bicharacteristic variables are not independent, 
insofar as their very definitions, eqns. (9), imply that: 


C + D 


V 2 


E + F 


V 3 " 


G + H 


(12,a,b,c) 


6a=C-D = E- F = G- H 


(12,d,e,f) 


so that 


F = -C+D + E 


(13) 


H = - C + D + G. (14) 

Therefore, any numerical solution obtained by integrating eqns. (8), directly, 
would lead to a "nonuniqueness" in the value of the speed of sound a. 
Furthermore, the right-hand sides of eqns. (8), namely the y^ coefficients, 
contain spatial derivatives of the velocity components, which are not 
associated with the convection of physical disturbances and are therefore 
likely to reduce the accuracy of the spatial discretization, if not the 
stability of the integration process. For these reasons, as in [5,6] for the 
two-dimensional case, the following equivalent system is obtained by taking a 
complete set of appropriate linear combinations of eqns (8): 
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3 C 3D 3 

c t + D t + <v i +a) js; + (v r a) + v 2 ss; (c + D) 


+ v 3 75; (C + D > ‘ - 2o i + 2e i 


( 15 ) 


E t + F t + v i Is; (E + r> + <v 2 +a > H; + H; 


+ V 3 % < E + F > ‘ -2ct 2 + 2S 2 


(16) 


3 3 3 n 

G t + H t + V 1 8^ (G + H) + V 2 (G + H) + (v 3 +a) 3^ 


, f . 9H 
+ (v 3~ a) 3Z 


- 2a 3 + 2g 3 


(17) 


¥ c t - \ * E t - F t + G t - H tl + < v i^ ) Hr - < v " a > 


3D 

i 9s i 


i * , \ 3E f v 3F , f . \ 3 G f N 3H 0 

+ (V 2'* a) 55; " <V 2‘ a) 55; + (v 3' te) 35; - (v 3‘ a) 55; ’ - 2C 


( 18 ) 


C t - r t - E t + F t ‘ 0 


(19) 


C t ' D t " G t + H t " °* 


( 20 ) 


It is noteworthy that eqns. (15-17) are simply the three components of 
the momentum eqn. (2), expressed as the sums of two compatibility conditions 
of two bicharacteristic variables along their bicharacteristic lines, whereas 
eqns. (18-20) all coincide with the continuity eqn. (1). Also, eqns. (19) and 
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(20), after a straightforward integration with respect to time, identically 
reproduce eqns. (13) and (14), so that they effectively reduce the number of 
dependent variables from six to four and any numerical integration of eqns. 
(13-18) will guarantee a unique solution for the physical variable a. 
Finally, the right-hand sides of eqns. (15-18) are seen to contain only 
source-like terms which do not involve spatial derivatives of the dependent 
variables and, therefore, are not likely to deteriorate the accuracy of any 
numerical method using windward differences for the "total " derivatives of the 
bicharacteristics variables. 

For these reasons eqns. (13—18) are considered the "most appropriate" 
three-dimensional lambda-formulation equations for a general orthogonal 
coordinate system and will be the basis for all of the numerical methods 
proposed in this study. 


NUMERICAL METHODS 

The governing equations (13-18) are discretized and linearized in time 
using the delta form [11,5,6] to give 

+ (u+a) n AC + (u-a) n AD + v n (AC 4- AD 1 + w n (AC 4- AD ) 

At At x x v y y J K z z J 

= . (u+a) n C* - (u-a) n - v n (C y + B y ) n - w n (C z + D z ) n (21) 

If- + + uI, ( AE x + AF X ) + ( v+a ) n AE y + (v-a) n AFy + w n (AE z + AF z ) 


= -u n (E + F ) n - (v+a) n E n - (v-a) n F n - w n (E + F ) n 
x x y y z z 


( 22 ) 
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^ H + u n (AG + AH ) + vYaG + AH ) + (w+a) n AG + (w-a) n AH 
At At x x / ^ y y z z 


= -u n (G v + Hj n - v n ( G v + H J n - (w+a) n g" - (w-a)“ H 


( 23 ) 


i / *£ 

3 Ut 


AD AE AF AG AHi , . , .n . .n . %n 

— -r-r + -r-r ~ -r-r + -r-r ~ T-rj + (u+a) AC - (u-a) AD + (v+a) AE 
At At At At At J x x \ 


(v-a) n AF^ + (w+a) n AG^ - (w-a) n AH^ = -(u+a)° + (u-a) n d” 


- (v+a) n E^ + (v-a) n f" - (w+a) n g” + (w-a) n H* 


AF = -AC + AD + AE 


(24) 


(25) 


AH = -AC + AF + AG, (26) 

where Cartesian coordinates have been used for simplicity, At is the time 

n 

step, AC = C - C (the superscripts n + 1 and n indicating the new and old 
time levels t n+ l _ t n + At and t n ) etc. Equations (21-26) constitute a 

first-order-accurate implicit time discretization of the corresponding 
differential problem; eqns. (21-24) are then discretized in space, using 
windward differences to properly account for the direction of wave 
propagation, and the AF and AH unknowns are eliminated in favor of 
AC, AD, AE and AG by means of eqns. (25) and (26) to produce, together with 
appropriate boundary conditions, a large 4x4 block-sparse linear system of 
the type 


A f = b. 


(27) 
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For the case of a cubic integration domain having N gridpoints in every 
spatial direction, A is a square matrix of order having only seven 

nonzero diagonals of 4x4-matrix-elements , f is the unknown vector having 
four-element-vector-components and b is the known coefficient vector. It is 
noteworthy that in previous works [5,6] a second-order-accurate time 
linearization was employed. However, due to the use of a backward Euler time 
discretization, eqn. (27) is only first-order-accurate in time anyway. 
Moreover, the present linearization, coupled with windward difference 
approximations for the left-hand sides of eqns. (21-24), leads to a diagonally 
dominant matrix A and has been verified to increase the stability of all the 
implicit methods later proposed in this study. It is also noteworthy that 
second-order-accurate, three-points windward differences can be used to 
approximate the right-hand sides of eqns. (21-24) so that, if the flow reaches 
a steady state, the final solution is second-order-accurate [5,6]. 

The main reason to employ an implicit method is to remove the CFL 
stability restriction, thus improving the efficiency of the calculations. 
However, a direct solution of problem (27), even if feasible, is certainly 
impractical. Therefore, the matrix A will be replaced by a matrix B which 
is easily invertible and is a first-order-accurate approximation (in time) 
of A. 

A Block-Explicit Method 

The simplest first-order-accurate approximation to A can be obtained by 
dropping all but the time-derivative terms in the left-hand sides of eqns. 
(21-24). The resulting matrix B is diagonal and a simple 4x4 linear 
system needs to be solved at every gridpoint to provide the local AC, AD, 

AE, and AG values. Furthermore, eqns. (21-24) can be rearranged to give 
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2AC 

At 


AC AD 
At At 


AC , AD 2AE 
At At At 


AC AD 2AG 

At At At 


= RHS(21 ) + RHS(24) 


= RHS(21) 
= RHS(22) 
= RHS(23) 


(28) 

(29) 

(30) 

(31) 


(where RHS(21) is a shorthand notation for the right-hand side of eqn. (21), 
etc.) so that every element of B is a lower triangular matrix which can he 
inverted directly. The present BE method has been developed mainly for 
assessing the efficiency of various implicit methods; however, due to its 
extreme coding simplicity, it could very well be a useful tool by itself, 
especially if implemented on a vector computer. 


A Block -Alternating-Direction-Implicit (BADI) Method 

An ADI technique has been developed, which is the direct extension to 

three-dimensional problems of the method of Refs. [5] and [6]. A three-sweep 

ADI process is used to solve problem (27) approximately. At the first sweep 

the t and x derivatives in the left-hand sides of eqns. (21-24) are 
evaluated implicitly, whereas the y and z derivatives are evaluated 
explicitly. At the second and third sweeps the t and y and the t and 

z derivatives are evaluated implicitly so that A is approximated by the 
product of three 4x4 block-tridiagonal matrices. In practice, at every 
sweep of the BADI method a 4x4 block-tridiagonal system of order N has to 

o 

be solved along each line of the computational grid, so that 3N such 

systems need to he solved at every time step (i.e., to solve eqn. (27) 
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approximately) . With respect to two-dimensional flow problems, the present 
ADI method is less competitive as compared to a standard explicit method, for 
two reasons; the block size of the tridiagonal systems increases from three 
to four and, more importantly, the number of tridiagonal systems to be solved 
at every time step grows from 2N to 3N 2 . Actually, for the simple problems 
later considered in this study the computer time per step for an ll 2 mesh 
was found to be about 30 times greater than that required by the BE method. 
More efficient implicit methods need therefore to be devised for the three- 
dimensional lambda-formulation equations. 

A Block-Line-Oauss-Seldel (BLGS) Method 

Classical relaxation methods have been recently employed with 
considerable success in connection with "upwind schemes" for the one- and two- 
dimensional Euler equations [7,8,12]. Here an obvious choice, leading to a 
reduction of the computer time per step to about one third, is to employ a 
single step 4x4 block— line— Gauss— Seidel method; all of the time and x 
derivatives in the left-hand sides of eqns. (21—24) are evaluated implicitly 
together with the diagonal contributions of the y and z derivatives, so 
that only N 4x4 block-tridiagonal systems (of order N) have to he solved 
at every time level. By accounting for the previously evaluated 
nontridiagonal entries explicitly, the matrix A is effectively replaced by 
its three main diagonals plus its two additional nonzero lower diagonals. 
Furthermore, the ordering of the solution process is changed at every time 
step so as to account for the two additional nonzero upper or lower diagonals, 
alternately. 
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A Block-Point-Gauss-Seidel (BPGS) Method 

By taking to its extreme the logic behind the previous method, an obvious 
choice presents itself; that is, to replace the matrix A with its lower or 
upper triangular part. In eqns. (21—24) the diagonal contributions are 
accounted for implicitly and the previously evaluated off-diagonal 
contributions are brought to the right-hand sides of the equations and 
accounted for explicitly. At every gridpoint location a 4x4 linear system 
needs to be solved as in the BE method; however, due to its variable 
coefficients, the local 4x4 matrix cannot be triangularized and a complete 
Gauss-Jordan elimination procedure, using diagonal pivot strategy, has been 
employed (here as well as to solve the local linear systems within a general 
block-tridiagonal inversion routine in all of the present implicit methods). 

A Simplified— Line— Gauss— Seidel (SLGS) Method 

From their very definitions (eqns. (9)) as well as from their compati- 
bility conditions (eqns. (8)) it appears that the waves associated with the 
bicharacteristic variables C and D mainly propagate in the x direction, 
whereas the E and F waves and the G and H waves mainly propagate in 
the y and z directions, respectively. Therefore, it would seem 
appropriate to devise a numerical method exploiting such a property of the 
compatibility eqns. (8), as done in [7,8] for the case of one- and two- 
dimensional flows. However, whereas Moretti [7,8] integrates the compatibilty 
conditions directly, here eqns. (13-18) are preferred for the two reasons 
previously discussed. In conclusion, the following simplified line-Gauss- 
Seidel method is proposed here: Equations (21) and (24) are solved coupled 
together for the AC and AD variables by means of a line-Gauss-Seidel method, 
implicit in the x direction, so that a 2x2 block-tridiagonal system of 
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order N has to be solved at every y • and z k grldpoint location. 
Equations (22) and (23) are then solved by means of line-Gauss-Seidel methods 
implicit in the y and z direction, respectively, so that 2N 2 additional 
scalar tridiagonal systems need to be solved. Obviously, equations (25) and 
(26) are used to eliminate AF and AH from eqns. (21-24) and all of the 

terms already evaluated at any level of the computation process are 
accounted for in the right-hand sides of the equations. Furthermore, since 
the pressure eqn. (24) does not have a main direction of propagation, it is 
coupled to eqn. (22), to evaluate AE and AF implicitly in the y direction, 
and to eqn. (23), to evaluate AG and AH implicitly in the z direction, 
at successive time steps. 

It is noteworthy that, in general, the matrix A contains all the boun- 
dary conditions, which are therefore accounted for with the level of impli- 
citness typical of every single method. However, for simplicity, in all of 
the present applications the exact solution of the continuum problem has been 
enforced at all boundaries to provide homogeneous boundary conditions for all 
of the incremental bicharacteristic variables. More general boundary 
conditions can be implemented as suggested in [6] and are not expected to 
cause any difficulty. 


RESULTS 

In order to test the proposed methods, a simple steady one-dimensional 
spherical source flow of air (y = 1.4) has been considered; for such a flow 
field the continuity and energy equations are given as 


14 


a 5 v r = c. (32) 

r 1 

0.2 v 2 + a 2 = c, (33) 

r 2 

v r being the radial velocity component and r the radial distance from the 
origin. All of the calculations have been performed using a Cartesian 
coordinate system inside the unit cube such that: 2 < x < 3; -.5 < y < .5; 

-.5 < z < .5. Three flow conditions have been considered: the subsonic flow 

corresponding to c^ = 3.2 and C£ = 1.128 and the supersonic and transonic 
flows corresponding to c^ = 4.2 and C£ = 1.2205. In the last case, an 
isentropic shock at r = 2.15 separates a supersonic region (for r < 2.15) 
from a subsonic one (for r > 2.15). The exact solution for the 
bicharacteristic variables has been imposed at all boundaries (the six sides 
of the computational cube) and a flow field having the exact values for u 
and a and zero v and w has been used as a suitable initial condition. 
The solution was advanced in time by means of any of the proposed methods 
using a constant (in time) and uniform (in space) value of At, until the 
average absolute value of AC at all interior points was less than 10 . 

Due to the use of the delta approach, the final steady solution is the same 
for all of the methods. The computed Mach number distribution along the x 
axis is plotted in Fig. 1 for the three flow cases versus the exact solution, 
for Ax = Ay = Az = .1. 

The solution is fairly good for the subsonic and supersonic case and 
qualitatively correct for the transonic one. In particular, the shock is 
captured in the correct mesh interval and no wiggles are present in spite of 
the absence of any additional dissipation. However, for shocks as strong as 
that given in Fig. 1, a shock-fitting procedure is warranted. 
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The computations were performed on a CDC Cyber 175 computer using two— 
windward differences for all of the spatial derivatives. 

The main purpose of this paper was to devise "efficient” implicit methods 
for the three-dimensional lambda-formulation equations. Therefore the 
performance of all of the present methods are given in Table I as the values 
of the At leading to the fastest convergence and the corresponding number of 
time steps (K) and CPU seconds. 


TABLE I 



Subsonic 

Flow 

Supersonic 

Flow 

Transonic 

Flow 

Method 

At 

K 

CPU 

At 

K 

CPU 

At 

K 

CPU 

BE 

.03 

178 

57 

.03 

62 

21 

.03 

403 

128 

BAD I 

.3 

35 

361 

.3 

38 

398 

.3 

86 

887 

BLGS 

>5 

18 

75 

>2 

11 

46 

>10 

71 

290 

BPGS 

>5 

29 

49 

>2 

11 

20 

>10 

76 

126 

SLGS 

2. 

22 

34 

.4 

27 

41 

>10 

99 

147 


From Table I the following conclusions can be drawn. For the supersonic 
flow case the BE and BPGS methods are clearly superior; this is obvious 
insofar as there is no upstream propagation in the x direction and thus the 
implicit methods use most of the CPU time accounting for zero entries. In 
terms of the number of iterations, the performance of the BLGS and BPGS 
methods are identical as they should be (all of the x derivatives being 
approximated with backward differences). For the more relevant transonic and 
subsonic flow cases the BLGS method always requires the smallest number of 
iterations to converge; however, the BPGS, SLGS and BE methods are the most 
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efficient ones, whereas the BADI method is consistently the least competitive 
one. It is noteworthy that all of the Gauss-Seidel methods are very robust 
insofar as they maintain a quasi-op timal convergence rate over a wide range 
of At values. Among the three "best methods," the BPGS and the BE methods 
are considerably simpler to code and require less computer memory, a very 
critical resource when dealing with three-dimensional problems. Therefore, 
preliminary studies have been conducted to assess the influence on their 
convergence rate of the mesh size and of second-order-accurate discretization 
for the nonincremental terms of the governing equations. The two methods 
converge in a number of iterations which is roughly inversely proportional to 
the step size (e.g., for the subsonic flow problem convergence is reached 
after 261 and 52 iterations for a 17 3 mesh and after 309 and 59 iterations for 
a 21 3 mesh, for the BE and the BPGS methods, respectively). However, it is 
noteworthy that for these calculations (performed on a VAX 11/750 computer) 
the BE method required about 2.5 more CPU time than the BPGS method. This 
indicates that the solution routine for the local 4x4 linear systems used in 
this study works less efficiently on the Cyber computer than the one used on 
the VAX and that the superiority of the BPGS method over the BE one is 
potentially greater than it actually appears from Table I. Also, the use of 
second-order-accurate differencing seems to deteriorate the convergence rate 
of the BPGS method less than that of the BE method. Finally, the superiority 
of the BPGS method (with respect to the BE method) is expected to increase 
even further by using a variable At [5,6,12] and when more general boundary 
conditions are employed; this, because the additional work will be relatively 
greater for the simpler BE method. 
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In conclusion, the BPGS method appears to be the most promising technique 
for solving three-dimensional compressible flow problems, by itself, or as a 
robust smoother within a more general multigrid procedure. However, both the 
BLGS and SLGS methods proposed in this study appear to be very promising 
alternatives to the ADI method of Refs. 5, 6 for solving two-dimensional 
steady flows, for which they are likely to outperform even the present BPGS 
method. 
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Figure 1. 


Numerical (symbols) versus exact (solid lines) solutions for 
spherical source-flows. 
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